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A hybrid method is developed based on the spectral and finite-difference methods 
for solving the inhomogeneous Zerilli equation in time-domain. The developed hybrid 
method decomposes the domain into the spectral and finite-difference domains. The sin- 
gular source term is located in the spectral domain while the solution in the region 
without the singular term is approximated by the higher-order finite-difference method. 

The spectral domain is also split into multi-domains and the finite-difference domain 
is placed as the boundary domain. Due to the global nature of the spectral method, a 
multi-domain method composed of the spectral domains only does not yield the proper 
power-law decay unless the range of the computational domain is large. The finite- 
difference domain helps reduce boundary effects due to the truncation of the computa- 
tional domain. The multi-domain approach with the finite-difference boundary domain 
method reduces the computational costs significantly and also yields the proper power- 
law decay. 

Stable and accurate interface conditions between the finite-difference and spectral 
domains and the spectral and spectral domains are derived. For the singular source term, 
we use both the Gaussian model with various values of full width at half maximum and 
a localized discrete (5-function. The discrete (5-function was generalized to adopt the 
Gauss-Lobatto collocation points of the spectral domain. 

The gravitational waveforms are measured. Numerical results show that the devel- 
oped hybrid method accurately yields the quasi-normal modes and the power-law decay 
profile. The numerical results also show that the power-law decay profile is less sensitive 
to the shape of the regularized (5-function for the Gaussian model than expected. The 
Gaussian model also yields better results than the localized discrete (5-function. 

Keywords: Dirac (5-function; finite-difference method; Spectral method; Hybrid method; 
Head-on collision of black holes ; Zerilli equation. 
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1. Introduction 

Since the classic works by Frank Zerilli in early 70 's on the particle falling 
in a Schwarzschild geomet ry, a lot of research and study has been performed on 
this fundamental problem I 4 | 5 | 6 | 7 | 8 | 9 | 10 | ll | 1 2j One of the earliest computational 
calculations was made by Press and his co-workers, which is now known as DRPP 
calculation ^ on the radiation emitted when a particle starting from rest at infinity 
falls into a non-spinning black hole. The collision of two black holes is, in principle, 
one of the most efficient mechanisms for the generation of the gravitational waves. 
In recent years the extreme mass ratio limit of the binary system has been a special 
focus of research in gravitational physics. Extreme-Mass-Ratio Inspirals (EMRIs) 
are one of the main sources of the gravitational waves for t he gravitational wave 
detectors, such as Laser Interferometer Space Antenna (LISA)^!^! EMRIs are binary 
systems composed of a stellar compact object (SCO) with a mass, m in the range of 
in = 1 ^ 10^ Mq inspiralling into a massive black hole (MBH) with a mass, AI in 
the range of M = 10* ~ IO^Mq located at the galactic center. Thus, the mass ratios 
involved are /i := ra/M ^ 10"'' — 10~^. During the slow inspiral phase the system 
is driven by the emission of gravitational radiation, the general features of which 
are now well understood. Press showed that there is always an intermediate stage 
where the ringdown is dominated by a set of oscillating and exponentially decaying 
solutions, quasinormal modes (QNMs) whose spectrum depends only on the mass 
of the black hole and the multipole- moment index I of the initial perturbation 
This regime is followed by a power-law tail decay due to backscattering. 

For the EMRI, the small companion black hole is modeled as a point particle, and 
the problem can be framed by using the black hole perturbation theory. Moreover, 
as the first approximation, the point particle follows the geodesies in the space- 
time of the central black hole. The frequency-domain approach to this problem has 
achieved many remarkable results. Specifically the accurate (< 10~*) determina- 
tion of the energy flux of gravitational waves was obtained in the frequency-domain 
11^. However, the frequency-domain approach can take long computational time and 
lose accuracy for non-periodic orbits (for example, parabolic orbits, orbits with high 
eccentricity or decaying orbits). The time-domain approach seems better suited for 
such orbits^. For the time-domain approach, the finite-difference (FD) method is 
one of the most popular numerical methods. The FD time-d om ain methods, how- 
ever, suffer from the relatively poor accuracy at the moment unless a very high 
computational resolution is used. The main reason is the point particle approxi- 
mation, i.e. the approximation of the singular source terms. Various approaches to 
this issue have been attempted, including the regularizing the Dirac (5-function us- 

ing a narr ow Gaussian distribution and also using more advanced discrete 5-model 
171161191201 

Anoth er approach of the EMRI problem is to use the spectral (SP) method 

| 7 | 17 | 18 | 9 | g^j. ppgvious work, we used the spectral method to solve the inhomo- 
geneous Zerilli equation in time-domain and obtained good results. But the proper 
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power-law decay was not observed In early time the solution agrees with the 
established solution but in very late time the solution is contaminated by the small- 
scale oscillations. These oscillations are likely due to the artificial truncation of the 
computational domain. 

In this work, we continue our previous research with the spectral method in 
order to obtain the proper power-law decay. For this, we developed the multi-domain 
hybrid method. The multi-domain method hybridizes the spectral method and the 
high-order finite-difference method. The spectral domain is also split into many 
sub-domains, each of which is also a spectral domain. The main advantage of the 
multi-domain method is that the computational costs can be significantly reduced 
by reducing the order of the interpolating polynomial in each sub-domain and the 
parallelization becomes robust. A fundamental reason for considering the multi- 
domain method is also to reduce the boundary effects due to the artificial truncation 
of the computational domain for obtaining the proper late time decay profile of the 
gravitational waveforms. In order to obtain the proper power-law decay, the outer 
boundary needs to be placed afar, in general. However, having a large size of the 
computational domain increases the computational costs significantly. In this work, 
we add the finite-difference domain as the boundary domain. The spectral method 
is a global method and it is highly sensitive to the boundary effects. To prevent 
the "fast" propagation of these boundary effects, we use a local method instead 
as the boundary domain, such as the finite-difference domain. By doing this, we 
obtain the proper power-law decay while having the computational costs reduced 
and also exploiting the accuracy of the spectral method. To patch each sub-domain 
with others, we derive the accurate and stable patching conditions. For the spectral 
and finite-difference sub-domains, we show that the resolution across the interface 
needs to be closely uniform. Otherwise, the CFL condition becomes strict. 

For the singular source term, we use both the Gaussian (5-function method and 
the discrete (5-function method. For the Gaussian method, we change the shape 
of the Gaussian profile to mimic the (5-function. For the discrete (5-function, we 
generalize the discrete (5-function developed by Sundararajan et al. '^^ into the 
one on the non-uniform grid. 

We provide numerical results that show the efficiency and robustness of the 
proposed hybrid method. Using the hybrid method we could obtain the proper 
power-law decay with the Gaussian approximation model. We use various shapes of 
the Gaussian profile and found that the result is insensitive to the shape. That is, 
even a broad profile, which results in a smooth solution, yields the power-law decay 
successfully. With the smooth solution, the spectral method does not need to use 
the filter operation, which increases the computational efficiency further. We also 
obtain the power-law decay with the discrete (5-function model, but the computed 
slope was not accurate, which may imply that the discrete (5-function model yields 
correct results only on uniform grids. 

This paper is organized as follows. In Section 2, we briefly describe the finite- 
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difference and spectral methods. For the finite-difference method, we used the 4th- 
order method. For the spectral method we use the Chebyshev spectral collocation 
method based on the Gauss-Lobatto collocation points. In Section 3 we describe 
the discrete (5-function on non-uniform grids. Section 4 explains the Zerilli equation 
briefly. In Section 5, we describe the proposed hybrid method in details. We derive 
the stable and accurate interface patching conditions. Boundary conditions arc de- 
scribed in Section 6. In Section 7, we discuss the stability of the hybrid method. In 
Section 8, numerical results are provided. In Section 9, a brief summary and future 
work are explained. 



2. finite-difference cind spectral methods 
2.1. finite- difference method 

In this work we consider both the 2nd and 4th-ordcr finite-difference method. The 
2nd-order finite-difference method for the spatial derivatives are well known and 
we omit those formulae. Instead we briefly explain the 4th-order finite-difference 
method. For the 4th-order method we will derive the formula when the grid is 
non-uniform. This is because, in the spectral domain, we use the Gauss-Lobato col- 
location points, which are not evenly distributed and we need the finite-difference 
formula for the boundary conditions in the spectral domain. Also we define the mod- 
ified fiux at the SP-SP interface, which also requires the finite-difference formulae. 
Details of these are described in section 5. 

4th- Order finite- difference method: Uniform grids. Let h be the grid spacing in 
the finite-difference domain and let ^ be a; < ^ < a; -|- /i. The standard Ath-ovAev 
derivatives are given by 

Centered difference: 

Off-centered (1 - point) difference: 
Off-centered (2 - points) difference: 

At the left x = x~ and right x = a;+ boundaries, we used the 2nd-order difference 
method. 

X = x~: x~ < ^ < x~ -\-h 

u'{x-) = -M^-)+M^-+h)-uix-+2h) ^ 2^.„(3)(^)^ 

where x~ < < x~ + h. 



Similarly 
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X = X~^: — h < < 

- 4u(a;+ - h) 
2h 

where x~^ — h < ^ < x~^. 



U'{X+) = Mx^)-^u{x^-h)+u{x+-2h) ^ 2^2„(3)(^)^ (5) 



The centered 2nd-order derivative is given by 

u"{x) = + 16n,-i - 30n, + IQu,^, - u,^, _ _L/,4„(6)(^). (g) 



2.2. Spectral method 

For the spectral method, we use the Chebyshev spectral collocation method based on 
the Gauss-Lobatto collocation points. The Chebyshev spectral collocation method 
seeks a solution in the Chebyshev polynomial space by the Chebyshev polynomials 
Ti{x) as 

N 

UN{x,t) = PNU{x,t) = ^Ui{t)Ti{x), 
1=0 

where Ti{x) is the Chebyshev polynomial of degree I and w; the corresponding 
expansion coefficient. The commonly used collocation points are the Gauss-Lobatto 
quadrature points Xj given by 

Xj=cos{^j), yj = o,---,N. 

These collocation points belong to [— 1 1] and so we need the proper mapping to fit 
in our computational domain. The required mapping is 

h — a , , , 

where [a, h] is the original computational domain. The expansion coefficients are 
given by 

2 ^ 1 

where c„ = 2 if n = , and c„ = 1 otherwise. 



2.3. Spectral filtering method 

We also use the spectral filtering method to minimize the possible non-physical high 
frequency modes. The oscillations with the spectral method possibly found near the 
local jiimp discontinuity and also generated due to inconsistent initial conditions 
propagate through the whole domain. Our filtered approximation is given by 

JV 

ul{x,t)=Y,<3/N)ujTj{x), 

3=0 
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where a{l/N) based on the exponential fiher is the fiher function according to^I^H^. 
Our filter matrix S is given by 

1 ^ 2 

where p is the order of filteration and e is a constant. The filtered solution at x = Xi 
is given by 

AT 



U 



N 



[xi) = ^ SijUj, Vi = 0, • • • ,N. 



3. Approximation of singular source 

3.1. Discrete S -function with non-uniform grid 

Discrete (5-function on a uniform grid has been derived in In SP-FD approach, 
the singular source term is always located ins ide the spectral domain. Since in the 
spectral domain the grid is non-uniform ^i^, we need to redefine the discrete d- 
function on the non-uniform grid. (5-function which exists aX x — a, Xk+i < a < 
Xk+2, by using the following relation 



^ hif{xi)Si = /(a) 

i 

{a - Xk+i){a - Xk+2){a ~ Xk+3) 



(Xfe - Xk+l){Xk - Xk+2){xk - Xk+z) 

[a - Xk){a ~ Xk+2){a - Xk+3) 



{Xk+l - 


Xk)ixk+1 


- Xk+2){xk+l - Xk+s) 


[a 


- Xk)ia - 


Xk+i)ia - Xk+3) 


{xk+2 - 


Xk)ixk+2 


- Xk+l){xk+2 - Xk+3) 


(a 


- Xk){a - 


Xk+i)ia - Xk+2) 


iXk+3 - 


Xk)iXk+3 


- Xk+l)iXk+3 - Xk+2) 



fixk) 

f{xk+i) 

fixk+2) 

f{xk+3)- 



Equating the coefficients of f{xi) from both sides yields 

' {a-Xk + i)ia-Xk+2)ia-Xk+3) ^ 

{xk + l-Xky'{xk + 2-Xk)(xk + 3-Xk) 

(a-Xk)(a-Xk+2)(a-Xk+3) „i „ 
, , ,1 Xk+l 



X. — ) {xk+l-Xk)ixk + 2-Xk + iy''(xk + 3-Xk + l) + (J) 

' 1 {a-Xk)(a-3:k+l)(o'-Xk+3) at Xi 

{Xk + 2-Xk)(xk + 2-Xk+l)(Xk + 3-Xk + 2y'' '^ + ^ 

{a-Xk){a-Xk+i)ia-Xk+2) ^ 

(xk + 3-^k)(xk + 3-Xk+l)ixk + 3-^k + 2)iXk + 4,-Xk + 3) "''+'^ ' 

For getting the first derivative of the S function we have 
^hJ{x^)^, = -/ (a) 

= - ^Kf ixi)Si 

STh X - f{xi-i) 

^ [Xt+i-Xi^l) 
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Equating the coefficients of f{xi) from both sides yields 

{xk + i~-a){o'~Xk+2}(a-'Xk+a) 



{a-Xk){a-Xk+2){a-Xk+3) 

iXk + l-Xk)^(Xk + 2-Xk){Xk + 2-Xk+l)(Xk + 3-Xk + l) 

{a — Xk + i){a — Xk+3) r a — Xk + 2 



{Xk + 2-Xk + l){Xk + 2~Xk) {Xk + l-Xk-l){Xk + l-Xk)(Xk + 3-Xk) 

a~Xk 1 

{Xk + 3-Xk + l){Xk + 2~Xk-i-l){Xk + 3-^k + 2}' 

{a-Xk){a-Xk+2) f a-Xk+1 

(Xk-\.3 — XkJr2)(Xk + 3 — Xk + l) (Xk + i~Xk + 2){x k + 3 — Xk)(Xk + 3-Xk + 2) 

a — XkJr3 ^1 

(Xk + 2-Xk){Xk+l-Xk)(Xk + 2-Xk + l)^ 

(o'-Xk)(a-Xk + i)(a-Xk+3) 

{xk + 3~-^k + l)(Xk + 2-Xk){Xk + 2-Xk + l)(Xk + 3-Xk + 2){xk + i-Xk + 3) 

_ (a-Xk)(a-Xk+l)(a-Xk+2) 

(Xk + 5-^k + i){Xk + i-Xk+2){xk + 3-Xk){Xk + 3-'>'k+l){Xk + 3-Xk + 2) 



at Xk-i 

clt Xj^ 

at Xk+i 

at Xk+2 

at Xk+d, 
at a;fc+4. 



(8) 



3.2. Gaussian S-function 

The Gaussian (5-function is defined by 

6{x) = ■ 

The shape of the (5-function depends on the fuh width at half maximum a. If a is 
small, S{x) is narrow and if cr has higher value then S{x) is broader. The value of /i 
depends on time and determines the position of the (5-function. 



4. Zerilli equation 

The lowest order perturbation theory of the initial Schwarzschild black hole space- 
time leads to the inhomogeneous Zerilli equation with even-parity Such an equa- 
tion describes the gravitational wave -0 in 1 -f 1 dimension given by the following 
second-order wave equation 

^tt=A^r'-Vi{r)^-Si{r,t), (9) 

where r* is the tortoise coordinate, Vi{r) the potential term and Si{r,t) the source 
term. For details see [3,10]. The tortoise coordinate, is given by 

r* = r + 2Mln{r/2M - 1), (10) 

where r is the physical coordinate and r > 2M. 

One can convert Eq (9) into a syste m (jf equations. In some previous works, Zerilli 
equation was solved in time-domain lEEl a,nd in all cases the authors converted 
the 2nd-order PDE to the system of Ist-order equations in space and time. We 
understand that some instabilities arose and to suppress those instabilities, they 
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introduced an auxiliary field which converted the equation to a coupled set of 1st- 
order equations in space and time. 

In this work wc do not convert the 2nd-order PDE to the system of equations. 
Instead we control the instability eflaciently by introducing new interface conditions 
between SP-SP interfaces and SP-FD interfaces. We also correct the numerical flux 
accordingly. 

5. Hybrid method 

The hybrid method is basically the domain decomposition method. Wc use 
the finite-difference and the spectral methods for the hybrid method. The one- 
dimensional single domain is decomposed into multi-domains. Each domain is car- 
ried by either the finite-difference method or the spectral method. 

Let n = [Re, Roo\ be the original domain. Suppose that is decomposed into M 
sub-domains and each sub-domain meets its adjacent sub-domain at its boundaries 
only such that 

M 

VL=\JVL\ Pi = d+VL\ or d-Q.\ j =i + l, 

i=l 

where is the ith sub-domain and 9+0 and d~fl denote the right and left bound- 
aries of the sub-domain. Let O* and fiy denote the ith sub-domain of the SP domain 
and the FD domain respectively. The ith sub-domain fi* should be either O* or Oy.. 
Let iV* -|- 1 be the total number of grid points of O' and let {xq, ■ ■ ■ , a;]^;} be the 
set of grid points in O*. The sub-domains are in order, indexed by the index i from 
1 to M. That is, i = 1 is corresponding to the leftmost sub-domain and i = M the 
rightmost sub-domain. In this work, by default, we let the last domain be the FD 
domain, i.e. = . The reason is to control the boundary effects arising from 
the outer boundary. We found some small oscillations arise from the outer boundary 
that contaminate the power-law decay when the multi-domain spectral method was 
used in late times. Since the FD method is the local method, we put a FD domain 
as the boundary domain to remove the boundary effect. 

The main development of the hybrid method is the stable patching algorithm at 
the domain interfaces between the SP and FD domains and the SP-SP domains. This 
involves the approximation of the derivatives across the non- uniform grid points, 
for which we consider the 4th-order FD method with the non-uniform grid and the 
polynomial interpolation method. 

5.1. 4th-order finite difference method: non-uniform grids 

The 4th-order centered FD method uses 5 points and the non-uniformity needs to be 
addressed in the FD approximation of the derivatives while the 2nd-order centered 
FD method does not need such non-uniform adaptation. The schematic diagram is 
given in Diagram 2. 
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2nd- order finite- difference derivative at the spectral domain boundary: FD do- 
main imposes the outflow boundary condition [dt -\- dx)u = Q B.t the right domain 
boundary. In the case that the last domain is a SP domain, we adopt the non- 
uniform FD method for the term of dx in the outflow boundary condition. 

Let hi = Xi — Xi-i and hi-i = Xi-i — Xi-2 and assume that hi and hi-i are aU 
non-zero constants and are not necessarily the same. With these definitions, let 

Hereafter let end and ghost denote the indices of the last grid points and the ghost 



grid points respectively. Then the derivative U'{x) at x = x^ ^ for fi*^ is given by 



X{X-l)hend • ^ ' 

4th- order finite- difference derivative at the interface across the adjacent spectral 
domains: At the domain interfaces across the SP domains, we use the 4th-order FD 
method for the 2nd-ordcr derivative (see Diagram 1). Let Xj—Xj-i = xj+i—Xj = hj 

, Xj — Xj-2 — Xjj^2 — Xj — hj+i and 

hj+i 



X = 



hj 



Then we have the 2nd-order derivative at x 

-Uj-2 + X^Uj-i - 2(A^ - l)Uj + X'^Uj+i - Uj+2 
/i2A2(A2 - 1) 

Using the above relation, the derivative a.i x = Xq for the SP domain Q.\ is given 

by 

-f/^-^ , + x'^uir^ 1 - 2(a4 - i)uk + x'^ui - m 

U {x^) /i2A2(A2 - 1) • ^^^> 



5.2. High-order interpolation 

For the interface patching between the FD domain and the SP domain, we use the 
high-order interpolation based on the Lagrange interpolation. For the interpolation 
at the ghost cells in the FD domain for the SP approximation in $7*, we use the 
4th-order polynomial. For example, for FD-SP patching we use the 5 points in the 
FD domain, These points are Zk = k = 0,1,2,3,4. The points for 

the ghost cells to seek are yi and 2/2 (two grid points marked with the symbol x in 
Diagram 2) and these are satisfying the following relations 

yi < 2/2 < x^'ti, 

X^iqi-l — y2 = x\ — Xq, 

J/2 - yi = 4 - x\. 
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We use the Lagrange interpolation to find the interpolation U{y) at y = yi, and 
y = y2 such that 

C/(y,) = ^C/-i(z,) n J-l>2. (14) 

For the interpolation at the ghost cehs in the SP domain for the FD domain (fi^, 
e.g.) we use the Chebyshev interpolation in the spectral domain ($7^+^). The grid 
points for the ghost cells in the SP domain are wi and W2 and these are satisfying 
the following relations (see Diagram 2) 

Xq^^ <Wi <W2, 

2+1 i i 

W2-W1 = X''j^i_^ - X^N--2- 

Since the FD domain has the uniform grids, we have wi — x^q^^ W2 ~ wi = 
x]^i — x^j^i_^. Similarly to the FD interpolation, the interpolation U (w) ai w — wi, 
and w = W2 is given by for j = 1 , 2 

U{w,) = U{xl+')Lk{i{w,)), J = 1, 2. (15) 

fc=0 

Note that ^ is the linear mapping from [oiq"^^, x^i+J to [—1 1]. Here Lk{w) is the 
Lagrange interpolation polynomial based on the Chebyshev polynomials given by 

^k(w) — 2 '■ 2 ' 

CfeiV*+i {w - Wk) 

where Cj = 2 if j 0, N^^^ and Cj — 1 otherwise and T^i+i {w) is the derivative 
of the Chebyshev polynomial of degree N^'^^ with respect to w. Since the polynomial 
order in each SP domain is low, the interpolation can be done quickly without using 
the fast transformation. 



5.3. Patching procedure 

To explain the patching procedure, we consider the following 2nd-order wave equa- 
tion, utt — Uxx — 0. Including the potential term is straightforward. 

SP-SP patching: Consider the two adjacent SP domains 51*'^^ and il* whose 
solutions are denoted by cj) and tp, respectively. (Hereafter, let and i/j denote the 
solution in the left and the right sub-domains, respectively.) The patching is based 
on Eq. at the domain boundaries only. For ri*+^ the left boundary value is 
updated based on the following 4th-order method: 
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where the superscripts j ' + 1, j and j — 1 indicate the time steps {j + 1)A<, 
jAt, and {j - l)At and ho = x\^'^ - x}^'^ and A = (a;^+^ - a;^+^)//io. Note that we 
let every SP domain has the same domain interval. Thus ho = x^j^i — x^^i_-^ and 

•^2 "^0 ~ "^N' •''JV'— 2" 

(§) (§) » '® — • • 



• Spectral grids, Q Spectral ghost cells 
Diagram 1: SP - SP patching. Grids at the SP and SP interface 

For Og the right boundary value is updated in the same way: 

,,+1 , ^^, -V>j,._, + AV^,_,-2(A^-l)^j,,+AV,-.^^, 
=2V'^,-V'^. +At /tgA2(A2 - 1) ■ ^ 



-e e o X ><§) « » A m — A«- 



O finite-difference grids, A finite-difference ghost cells 
• Spectral grids, x Spectral ghost cells 

Diagram 2: FD - SP patching. Grids at the FD and SP interface. 



FD-SP patching: For the FD and SP domain patching, wc use the polynomial 
interpolation for the ghost cells. For this case, the 4th-order method is used based 
on Eq. 6 for the uniform grid. Let tjj and (f> be the solutions in the FD domain fl^ 
and SP domain fl*"*"^, respectively. The value at the right boundary of fij is updated 
using the following: 



= 2V.;. - v^^-i + At'^^ — uh' — - — (18) 



7 
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where hf — x^j^i — x^j^i_-^ and are from Eq. ([15 



fe=0 



The value at the left boundary of ri*'^^ is updated in the similar way, but using 
the non- uniform formula Eq. (|16p . we have 



where ho and A are defined in the same way as in Eq. (jl6p . And the interpolation 
values are from Eq. ([H]) given by 

where = ^^N'^-k ^^'^ 

Vm = a;]v' ~ ^oA™, m = 0, 1. 

Flux correction : After the patching procedure for the FD-SP and SP-SP 
patching, we consider the flux correction. Since we set each SP domain to have 
the same resolution, it is obvious that 

cl>'{ghost) = ^[2), 

ip-' [ghost) = ^j-' (end — 1). 
We use the central (average) flux to get 

Flux = i[7/'^'+i(end) + </'^'+^(l)] 

= ]^[2ijj^ {end) - iP^-^{end) + 2<?!)^'(1) ~ (j}^-^{l)] 

+ i(^)2[20J(2) + 2ip\end - 1) - 2il^\end) ~ 2(/ij(l)] 

= ]^[2ij^{end) - f^-^{end) + 2<?!)^'(1) - <?!)^'"^(1)] + 

(A^)^ 0J(2)-0j(l) ^^{end) - il)^{end-l) 
h ^ h h ^ 

- -[2V^(end) - V^-i(end) + 20^(1) - ^^-^(l)] + ^-^[0.(1) - M^nd)]. 

If the solution is smooth, e.g. ip,(t) ^ then (t)x{^) — tpx(end). Then the flux 
reduces to 

Flux = i[2?A^'(end) - ^^-^(end) + 2(/)^'(l) - (jy^^^l)]. (20) 
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For example, at steady state ^^{end) = ip^~^{end) = ^^'^^{end) and so the flux is 
Flux= ^[tp^ {end) + (/>>{!)]. 

Similarly for the 4th-order method, 



X'^(l)^{end - 1) - 2{X^ - l)^ (end) + X'^^j^ {ghostl) - ip^ {ghost2)\ 



^+\end) = 2<P^{end) ~ <P'-\end) + T2T27T2 tt h-/-^ (end - 2) + 

/if A^(A^ — 1) 



and 



Af2 

^.+1(1) = 24>'\l) - + h2x2^x^_i) [-^i9host2) + 

X'^(jP{ghostl) - 2{X^ - 1W{1) + A^V'(2) - V-^^S)]. 

For SP-SP hybrid process, and at the jth interface {ghost2) = ^ {end — 2) , 
(j)j{ghosfA) = (f)i{end- 1) , i'i{ghosfA) = tpi {2) and il)^{ghost2) = ■^■'(S). We can 
write the central flux at the jth interface as follows 

pi^^^f^end) + r+^{l) 



2 

{2(f>i{end) - (t>>-'^{end) + 2i>i{l) - V'^-^(l)) 
2 



+ 



[-(jP{end-2) + 



2/i2A2(A2-l)' 
X'^(l^{end - 1) - 2(A'' - IW {end) + X^i!^ {ghostl) - 
ip^ {ghost2) + (j^{ghost2) + 

X^(IP{ghostl) - 2(A* - 1)V^(1) + AV'(2) - -0^(3)]. (21) 

If we use the central flux directly then 

ip^+'^{end) = 2'iP^{end) - ij^-^{end) + {At)^D'^^^, 
(/)^+i(l) = 20^(1) - 0^-^(1) + {M)'^D'^4P, 



(22) 



where D is the spectral differentiation matrix. 

Flux = ^[2tp^{end) - ilj^-'^{end) + 2(fP{l) - ^''^{1)] + 



{At) 



2 



[D'€nd + D'ct>i]. (23) 



If the function is not then D'^tpl^^^ ^ iP^^J Pqj. ^^^e FD-SP domain the resolu- 
tions arc usually different, i.c (fp {ghost) 7^ 0^ (2) and ip^ {ghost?) 7^ ih^ {end — 1). But 
we can expand (p^ and tp^ around the 2"^^ and {end — l)-th point respectively. 
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6. Boundary conditions 

For the boundary conditions at the boundaries 9f2 of the whole domain £7, we use 
the simple outflow condition based on the assumption that the potential term at 
dVl is negligible. That is, we use 

=0, r*ed-n\, 

+ =0, r*ed+nf. 

For the first-order derivative in the outflow equations, we use the 2nd-order finite- 
difference method, Eq. ^ for the FD domain $7^^ and Eq. (|lip for the SP domain 




,"+1 ,/,"+! 
N ' Vl 



7. Stability analysis 

O (s) 



O SPl, • SP2 

Diagram 3; SP - SP patching. Grids at the SP and SP interface. 



7.1. SP-SP stability analysis 

For stability analysis we consider SP-SP domain for example. This analysis is for 
two domains. Let us consider the two spectral domains denoted by SPl and SP2. 
Let 4> be the solution in SPl and ^/i be the solution in SP2 of the Zerrilli equation 
without the potential and singular source terms. Then we have 

= 2</>7 - + {MfDl^^l J - 1, • • • , iV - 1, 

= 2^^ - + {AtfDj^l J ^ 1, . . . , iv - 1. 

We assume that there is no boundary effect. Also we consider the equal length 
intervals and so Di = D2 — j^D where D is the original Chebyshev differential 
matrix and L is the length of the spectral sub-domains. Df is the sub-matrix of Df 
without the last and first rows and the first column and the sub-matrix of 

without last row, first row and first column. 
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Collecting all the above equations in matrix form yields, 







4>1 




(f'l 

J n— 1 
02 




= 2 


<t>l 

















{Df)N-l,N 




00 



1 _ 1 n 1 
1 1 



27F 









N 



2)N-1,N 





















rN-i. 



where N is the N—lx (N) null matrix. The previous matrix equation can be written 
in compact form 



= 2W" - W"--' + DW, 



n-l 



or 



W = I2nW". 



Eqs 23 and 24 can be written in matrix form 



Let 







2I2N + D —I2N 




W" 






. I2N . 








A = 


2/2Af + D —I2N 








_ I2N _ 





(24) 
(25) 



For the matrix stability, we have the spectral radius of A less than 1. i.e p{A) < 1. 



To check the stability we consider two spectral domains [—1,1] and [1,3]. For these 
two domains Di and D2 must be equal. Let hi be the minimum grid spacing in 
the left domain and similarly /12 for the right domain. Here we have hi = h2 = h. 
We found the following results for ^''(number of grid points in each domain) and 
At(time step) for which p{A) < 1. 
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Table 1. Comparison of number of grid points (A'^), time 
step (Af), minimum grid spacing (h). 



N 


At 


h 


16 


10-5 


7.810 xlO-3 


32 


10-6 


1.953 xlO-3 


64 


10-6 


4.883 xlO-4 


128 


10-7 


1.221 xlO-4 


256 


lo-'^ 


3.052 xlO-5 


512 


10-8 


7.629 xlO-^ 



7.2. SP-FD stabilty analysis 



O. 0.(S) 



hi h h h 
O SP, • FD 

Diagram 4; SP - FD patching. Grids at the SP and FD interface. 



At the SP-FD interface the grid resolutions between the adjacent sub-domains 
across the domain interface are different. So the grid distribution is non-uniform. 
The stable interface conditions derived for the uniform grid system are not enough 
and we need some conditions with which the spatial non-uniformity can be addressed 
properly [see 26 and references therein] . Let A'^i be the number of grid points in each 
SP sub-domain and be the number of grid points in the FD domain. Also assume 
Li be the length of each SP sub-domain and L2 be the length of the FD sub-domain. 
Then for stability we must have^^B 



2- .r- (26) 



Consider one spectral domain (SPl) and one FD-domain (FDl). Let (j> be the 
solution in SPl and ip be the solution in FDl of the simple wave equation without 
any potential term. Let number of grid points in SPl and FDl be iVi and N2 
respectively which satisfy eqn (25). We assume that there is no boundary effect. 
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Then from spectral collocation and 2nd order FD-method, we have 
<^^"+i = 24>^ - + (At)2i520^", J = 1, . . . , TVi - 1, 

V;;+i = 2^-," - + (A*)'^2V'i, J = 1, • • • , VV2 - 1. 

Where Di = j^D where D is the original Chcbyshcv matrix, D\ is the sub-matrix 
of D\ without the first and last rows and the last column and D2 is the 2nd-order 
FD-differentiation matrix, which is written by 

1-2 1 
1-2 1 



2/i2 



and D2 is the sub-matrix of D2 without the first and last rows and the first column. 
Collecting all the above equations in matrix form yields. 







4>i 








= 2 






in — l 













■l,JVi 








1 
1 

27? 



No 



-In 1 
n - 1 1 



JV2-1,JV2 














































Where Ni is the Ni — lx {N2) null matrix, and N2 is the N2 — lx (Ni) null matrix. 
The previous matrix equation can be written in the following compact form 



= 2W" - W""^ -h bW' 
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or 



VF"+^ = (2J + £>)M^" - 



(27) 
(28) 



Eqs 26 and 27 can be written in matrix form 







'2i + D -T 










i 







Let 

2/ + L> -/ 

Where I is the unit matrix of order Ni+N2. For the matrix stabiUty, we must have 
the spectral radius of A less than 1. i.e p{A) < 1. 

To check the stability we consider the SP-domain [—1 1] and FD-domain [1 2]. 
We choose Ni and N2 in such a way that they must satisiy Eq. (25). We found the 
following results for A''i (number of grid points in SP-domain), N2 (Number of grid 
points in FD-domain) and Af(time step) for which p{A) < 1. Some examples are 
given in the following table: 

Tabic 2. Comparison of number of grid points in SP do- 
main (Ni), number of grid points on FD domain {N2) and 
time step At. 



N2 



At 



8 64 10-5 

16 256 10-^ 

32 1024 10-''' 

48 2304 10-* 



The tabic implies that the more uniformity of the grid spacing across the SP and 
FD interface is achieved the larger value of time step could be used for stability. For 
our numerical experiments in the next section, we choose the geometric parameters 
for the sub-domains so that the grid resolution is uniform in the small neighborhood 
across the SP-FD interface. Here note that for the numerical experiments in this 
paper, we choose the geometric setting so that we have about At ^ 10~^. 

8. Numerical results 

For the numerical experiments, the computational domain is split into two main 

sub-domains, i.e. the SP and FD domains. The (5-function is located in the SP sub- 
domain for all time (see Diagram 5). The SP sub-domain is also divided into multiple 
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smaller SP sub-domains and each sub-domain communicates with its adjacent sub- 
domains by the interface condition described in Section 5. The last SP sub-domain 
communicates with the boundary FD sub-domain. 



Diagram 5: SP - SP .... - SP - FD subdomains. 



Since the spectral method is a global method, the boundary effects spread in- 
stantly throughout the domain. Thus, the SP method suffers from the boundary 
effects unless the range of the domain is large enough. However, the multiple SP 
sub-domains with the FD boundary domain help avoid the unphysical oscillations 
and reduce the computational cost. 

To determine the number of sub-domains and the order of the interpolating 
polynomial in each sub-domain, two aspects are considered. The first aspect is the 
resolution for the singular source term and the second one is the grid uniformity 
across the SP and FD domain interface. As explained in the previous section, the 
non-uniformity of the grid resolution near the SP and FD interface makes the CFL 
condition strict. Considering these aspects, we truncate the domain to [—300, 387.5] 
in r*. The SP sub-domain is [-300,20] and the FD sub-domain is [20,387.5]. The 
multi-domain SP method reduces the computational time significantly by reducing 
the size of the system but still achieves the desired accuracy. The typical setup in 
this work is that the number of sub-domains in the SP domain is {N) = 40 and each 
domain has the interpolating order of (n) = 48. The time-step is relatively large 
such as 10~^ and the Gaussian 5-function has a = 20. The waveform is collected 
at various values of r*. 

8.1. Gaussian S -function 

First we consider the case with the Gaussian model. The first result is in Figure 
1. The waveform is collected at r* = 137.6. The left figure in the top panel shows 

the ringdown profile which starts around the time t = 150. The right figiirc in the 
top panel shows the power-law decay. The two figures in the lower panel show the 



SP sub-domains 



FD sub-domain 
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two distinct phases. For the time approximately t — 150 ^ 300 the solution decays 
exponentially and it is oscillatory. This phase is the Quasi-Normal-Mode (QNM) 
ringing phase. After thi s ph ase the power-law decay starts. According to the seminal 
work by Richard Price the observer measures the late time perturbation field 
to drop off as an inverse power law of time, specifically as In the case of 
Schwarzschild black hole, n ~ 21 + 3, where I is the multipole moment of the 
perturbation field. In the context of our present work*-^, I = 2. The right figure in 
lower panel shows the power-law decay in logarithmic scale both in the horizontal 
and vertical axes. In this phase, the slope of the decay profile is about ^ The 
slopes of the power-law decay for various a and different reso lutions, are almost 
same and ~ —7. This agrees well with the theoretical result ^ i"^'""^. 

In Figure 2, we use the Gaussian J-function with a = 10. The waveform is 
collected at r* ~ 137.6. The top panel shows the ringdown profile and lower panel 
shows the QNM and the power-law decay. The similar QNM and the power-law 
decay are observed. Also we find that the same number of oscillations are observed 
in the QNM regime as in Figure 1. 

In Figure 3, we use the Gaussian model with a = 40. The waveform is collected 
at — 137.6. Note that the higher value of a is used for this figure and the Gaussian 
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200 300 400 500 2.3 2.4 2.5 2.6 2.7 

tM log, ftM) 



Fig. 2. The QMN and the power-law decay for Gaussian model with cr = 10. Rc = —300, 
iJoo = 387.5, ro = 1.5(2M). 

(5-function is much smoother than the previous two cases. We, however, found that 
the similar solution with the expected QNM profile and power-law decay profile was 
obtained. 

In Figure 4, we use the Gaussian S model with a ~ 50. The waveform was 
collected at — 137.6. With this high value of cr, the Gaussian (5-function is 
even smoother. This high level of smoothness of the singular source term made 
it unnecessary to use the spectral filter for the solution. For the previous figures, 
we used the spectral filter to regularize the solution due to the possible Gibbs 
oscillations. Without the filter operation, the computational time was reduced. This 
suggests that the desired decaying profile can be obtained by having the parameters 
in a clever way. For example, in Figure 5, we show the polynomial order n that we 
used for figures with different values of a to obtain the desired QNM and the power- 
law decay. The parameter values in the figure are (n, a) = (96, 10), (48, 20), (32, 40), 
(24, 50). In our current research we didn't attempt, but it would be an interesting 
study to investigate the optimal configuration of these parameters. 

8.2. Discrete 5-function 

We repeated the above experiment with the discrete i5-function. By definition, the 
discrete 5-function is localized and its shape changes with time because the spectral 
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t/M log, ftM) 



Fig. 3. The QMN and the power-law decay for the Gaussian model with a = 40. Re = —300, 
= 387.5, ro = 1.5(2A/). 

grid spacing is non-uniform. 

The QNM and the power-law decay profiles are depicted in Fig. 6. But the slope 
of the power-law tail does not match with that for the Gaussian model. For example, 
in Fig. 6 we considered the case that = —300, Too — 500 and the interface between 
the SP and FD domains is at = 150. We use 40 SP sub-domains with n = 48. The 
waveform was collected at — 250. The decay rate is slower than expected. The 
estimated slope is ~ —2 while the slope with the Gaussian model was close to —7. 
Although we increase the resolution, n, no significant improvement was observed. 
Fig. 7 shows the QNM and the power-law tail with the discrete (5-function. Here we 
used the very far outer boundary Too = 1970. The slope is still about —2. It seems 
that the discrete 5-function is good for the FD domain with the uniform grid but 
it is not ideal for the non-uniform grid with the SP domain. 

9. Conclusion 

In this work, the inhomogeneous Zerilli equation was solved numerically in time- 
domain, for which we developed a multi-domain hybrid method based on the Cheby- 
shev spectral collocation method and the 4th-order explicit finite-difference method. 
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Fig. 4. The QMN and the power law decay for the Gaussian model with a = 50. Re = —300, 
Roo = 387.5, ro = 1.5(2M). 



120r 
100 

80 

60 




20 

10 20 30 40 50 

a 

Fig. 5. The polynomial order n in each SP sub-domain with the value of cr for the Gaussian 
model. TV = 40. 

Using the developed method, we computed the waveforms for head on coUisions of 
black holes in one dimension. The in-falling black hole is modeled as a point source 
and consequently singular source terms appear in the governing equation. For the 
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Fig. 6. The QNM and the power-law decay for the discrete (5-function. Re = —300, Roo = 500, 
ro = 3(2M). 




X 10 



Fig. 7. The last few QNM modes and the long power-law tail for the discrete <5- function. 
Re = -30, iJoo = 1970, ro = 200(2A/). 



approximation of singular source terms, we used both the Gaussian (5-function and 
the discrete 5-function. For the stable and accurate approximation, we derived the 
interface conditions between the spectral and spectral domains and the spectral and 
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finite-difference domains. Tlie main approaclr introduced in tlris worfc is tire use of 
tlie finite-difference domain as ttie boundary domain. Witliout tire finite-difference 
domain as the boundary domain, the multi-domain composed of only the spectral 
sub-domains does not yield the proper power-law decay profile unless the range of 
the computational domain is very large. Using the multi-domain approach with the 
finite-difference boundary domain method, we could obtain the proper power-law 
decay profile with a relatively small computational cost. That is, the CFL condition 
is much relaxed and the location of the outer boundary of the computational do- 
main is not afar. Numerical results show that the hybrid method obtains a proper 
quasi-normal and power-law decay with the Gaussian 5-function approximation. 
Interestingly, even with a large value of tr, the proper power-law decay was ob- 
served. With the large value of a, the spectral filtering operation was not necessary, 
so the computational time was much reduced. The hybrid method with the dis- 
crete 5-function approximation, however, did not yield the proper power-law decay. 
The current study only considered the multi-domain spectral and finite-difference 
method with the (5-function residing in the spectral domain. We will investigate 
the optimal configuration of the multi-domain computational domain in our future 
research. 
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